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Force fluctuations in granular materials are investigated. A continuum equation is derived start- 
ing from a discrete model proposed in the literature. The influence of boundary conditions is 
investigated. For periodic boundary conditions the average weight is found to increase linearly with 
depth while it saturates to a constant value for absorbing boundary conditions, which models the 
existence of walls. The scale dependencies of the saturation weight, the saturation depth and the 
average squared fluctuations are obtained. The analytical results are compared with previous works 
and with numerical simulations in one dimension. 
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I. INTRODUCTION 

A wide variety of technical processes involve the stor- 
age and transport of granular materials. On the other 
hand, granular materials have very unusual properties 
which have intrigued researches in physics over recent 
years [0. Of great importance is the characterization 
of stress fluctuations in dry granular materials. Optical 
measurements in two and three dimensional Q arrays of 
granular materials have showed that the stress in packed 
granular materials is not distributed uniformly inside the 
medium, but is concentrated along " chains" . Further ex- 
periments dealing with force fluctuations in the bottom 
of relative small containers have shown an exponential 
distribution of vertical forces (weight) w 

Another interesting phenomena associated with weight 
fluctuations is the problem of arching. Arching refers to 
the observation that the average weight at the bottom of 
a column of grains saturates at a value W s as the depth t 
measured from the top of the column is increased, where 
W s shows large fluctuations when repeating this proce- 
dure with the same amount of grain. About one century 
ago Janssen [0 showed, using a simple argument, that 
W(t) = W s [l — exp(— t/t s )]. t s is the saturation depth 
which scales as t s ~ L, where L is the linear horizontal 
size of the container. 

More recently, Liu et al jj) proposed a simple discrete 
model (q-model) which describe the essential features of 
force fluctuations in granular materials. For instance, 
numerical simulations of the model in 2 + 1 dimensions 
give rise to an exponential distribution P(v) ~ exp(— Xv) 
for the normalized weight v — w/t, in agreement with 
experimental observations and numerical simulations of 
sphere packings (||| ■ Later Peralta-Fabi et al |tJ showed 
through numerical simulations in 1 + 1 dimensions and 
mean field (MF) analysis that the same model, but with 
different boundary conditions, explains the process of 
arching, obtaining the Janssen law but with W s ~ L 2 
and t s ~ L 2 . 

Generalizations of the q-model has been considered 
|p|-pT[. Some of them including a local slip condition 
which leads to anisotropies in the stress transmission M 



and introduces correlations ||. These correlations, nev- 
ertheless, do not alter the asymptotic exponential decay 
of the normalized weight distribution ||. On the other 
hand, other authors have introduced more realistic mod- 
els which takes into account the stress fluctuations in the 
horizontal cross-section ■ Some of their predictions 

are in agreement with those obtained for the (/-model, 
such as the exponential asymptotic decay of the normal- 
ized weight distribution characteristic of the g-model. 

We also count with some coarse-grained equations for 
the g-model QJl^l . Claudin et al jl2| obtained a diffusion 
equation for the weight with a multiplicative noise and 
derived some analytical results for the weight correlation 
function. However, their analysis was limited to pile con- 
figurations and silos with periodic boundary conditions. 
On the contrary, Peralta-Fabi et al the case when part 
of the weight is supported by the walls of the silo, which 
lead to a saturation of the average weight profile. Nev- 
ertheless, they neglected fluctuations in the local stress 
transmission. 

In the present work we focus our attention in obtaining 
a continuum model to describe weight fluctuations and 
the process of arching in granular materials, for both pe- 
riodic and wall boundary conditions. We take as starting 
point the q- model introduced by Liu et al jj| . We obtain 
a coarse-grained continuum equation of this model which 
is similar to the Edwards- Wilkinson equation, where the 
grain weight acts as an external force and noise is added 
to consider fluctuations in shape and orientation. Bound- 
ary and initial conditions are imposed following the def- 
inition of the discrete model. Two particular cases are 
analyzed, periodic and absorbing boundaries. In the case 
of periodic boundary conditions the average weight in- 
creases linearly with depth and is independent of the lat- 
tice size, while in the case of absorbing boundary condi- 
tions it saturates for t 3> t s ~ L 2 to W s ~ L 2 . Moreover, 
in both cases fluctuations around the average are of the 
order of L c , with C = (4 - d)/2. 

The paper is organized as follows. In section [I] we 
derive the coarse-grained equation for the (/-model and 
solve it for the case of periodic and absorbing bound- 
ary conditions. In section III our results are compared 
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with numerical simulations in one dimension, and with 
other works in the literature. Finally the summary and 
conclusions are given in section IV. 



II. THE MODEL 



B. Coarse-grained equation 



To derive a coarse-grained continuum equation of eq. 
(0) let us introduce the new random fraction Qij(t), such 
that 



A. Discrete model 

Liu et al || introduce a simple model (g-model) that 
describe some of the experimental observations in bead 
packs. The model assumes that the dominant physical 
mechanism leading to force fluctuations is the inhomo- 
geneity of the packing, which causes an unequal distri- 
bution of weight on the beads supporting a given grain. 
Spatial correlations in these fractions and variations in 
the coordination numbers of the grains are ignored. Only 
the vertical components of the forces are considered ex- 
plicitly. 

Consider a d + 1 dimensional array of grains in which 
each layer has L d beads. The total weight on a given 
bead is transmitted unevenly to N adjacent beads in the 
layer underneath. Specifically, a fraction of the total 
weight supported by the bead j in layer t is transmitted 
to bead i in layer t + 1. Thus, a site at depth t + 1 has 
weight Wi(t + 1), due to its own weight wq and to weights 
of neighbors at depth t, according to 



Wi(t+ 1) = w +y]qij(t)wj(t), 



(1) 



where the sum runs over the N neighbors, in layer t, of 
the site i, in layer t + 1. The fractions <?y(t) are taken 
to be random variables, independent except for the con- 
straint 



(2) 



where the sum runs over the N neighbors, in layer t+ 1, 
of site j, in layer t. 

To complete model definition, boundary and initial 
conditions must be specified. In the top layer (I = 0) 
the grains have no neighbors above and, therefore, they 
only support their own weight, i.e. u>i(0) — Wo- On the 
other hand, we can take different boundary conditions. 
We consider two cases, periodic and absorbing boundary 
conditions. Absorbing boundary conditions are more ap- 
propriate to describe the force fluctuations in silos where 
part of the weight is supported (" absorbed" ) by the wall 
containing the grains. The absorbing boundary condi- 
tions are implemented by simply imposing that boundary 
sites have zero weight. If a site near the boundary trans- 
mit part of its weight to a boundary site (i.e. to the wall) 
this weight fraction is "dissipated", while the boundary 
site never transmit weight to its neighbors because by 
definition its weight is zero. This boundary conditions 
arc different but equivalent to the one used by Pcralta- 
Fabi et al [Q. On the other hand, for periodic boundary 
conditions we recover the original g-model H . 



(3) 



which following (g) satisfies the constraint 

= 0. (4) 



Substituting in eq. (|I]) by eq. ||) and substrating 

Wi(t) to both sides it results that 

Wi(t + 1) - Wi(t) = — J^jwjjt) - Wi(t)} +WQ+ T]i(t), 



where 



%(*) = ^Qij(*K-(*)- 



(5) 



(6) 



rji (t) is a noise term associated with the random fractions 
Qij{t) and with weights at neighbor sites. In average the 
contribution of rji (t) should be zero because of the con- 
straint in eq. (^). 

Eq. (||) can be coarse-grained to obtain a continuum 
equation for the effective w(x, t). In the left hand side we 
have the discrete depth derivative, while the first term in 
the right hand side is the discrete Laplacian in the hori- 
zontal direction. After coarse-graining it result that 



T\7 2 w(x, t) + w + r)(x, t), 



(7) 



A and V are coarse-grained coeficients. A ~ a± and 
r ~ a?JN, where a± and <Z|| are characteristic lengths 
in the vertical and horizontal direction, respectively. 

In (1 + l)-dimensions eq. (|7) is actually very similar to 
the one derived in . Claudin et al [Q considered the 
explicit multiplicative nature of the noise. However, it is 
not clear how their analysis can be extended to (2 + 1)- 
dimensions, in particular the precise form of the coarse- 
grained noise is in this case not clear for us. Instead, we 
assume that T)(x, t) is a Gaussian noise with zero mean 
and uncorrelated in space, with noise correlator 



{ri(x, t)r]{x\t')) = 6(2 - x')A(t - t'). 



(8) 



A(t) is a monotonically decreasing even function. We do 
not know the precise form of A(t), but we assume it de- 
cays to zero beyond depth t c . Depending on the value 
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of t c we can obtain different behaviors. For t c ~ a± the 
noise will be also uncorrelated in the vertical direction, 
i.e. A(f) ~ 6(t). On the contrary, if t c is much larger 
than any characteristic depth then A(t) may be consid- 
ered constant. We expect the noise rj(x, t) to be strongly 
correlated in the vertical direction. According to eq. (g) 
the noise actually depends on weights at the different 
sizes, which should be strongly correlated from layer to 
layer. We have avoided the multiplicative nature of the 
noise but in compensation we must keep the correlations 
in the vertical direction. In any case we are going to con- 
sider different choices of A(t) and compare the results 
with numerical simulations. 



C. Average weight and fluctuations 

Eq. (j7j) can be interpreted as the equation of motion 
of an interface profile w(x, t), where depth plays the role 
of time, under an external force wo and annealed noise 
rj{x,t). This observation is very important because it 
shows that force fluctuations in granular media can be 
described through a more general framework, that if in- 
terface dynamics, which has been extensively studied in 
the literature In this context, eq. (^) is known as 
the Edwards- Wilkinson (EW) equation after 0. 

A central quantity of interest is the width Aw 
of the fluctuating "interface", given by Aw 2 (L,t) = 
(L d J Q d d [w(x,t) — (w)} 2 ), where (w) is the mean height 
of the interface (mean weight). In general the "surface 
roughness" Aw has the following asymptotic behavior 



Aw(L,t) = L^f(t/L z ), 



(9) 



where £ and z are the roughness and dynamic exponent, 
respectively. In the thermodynamic limit L — > oo, being 
linear, eq. (g) is readily solved via Fourier methods. Di- 
rect integration shows that, if the noise is spatially and 
temporally uncorrelated, the roughness exponent is given 
byC = (2-d)/2 H|. 

In this section we investigate the solution of eq. (Q) 
for different boundary conditions and noise correlators 
A(t). We are interested in the stationary solution and, 
therefore, the initial condition is irrelevant. 

Let first analyze the case of periodic boundary condi- 
tions. In this case the solution can be written as 



w(x,t) = -yt + y(x',t). 



(10) 



After substituting this expression in eq. (^) we obtain 
the following equation for y(x,t) 



\d t y{x,t) = TV 2 y{x,t) +r 1 {x,t). 



(11) 



The solution of this equation clearly has zero average, i.e. 
(y(x,t)) — 0, and therefore 



(w{x,t)) = -j-t. 



(12) 



The average weight thus increases linearly with t and no 
saturation is observed. 

On the contrary, for absorbing boundary conditions 
the solution cannot be proposed as in eq. (|l0|). In this 
case weight dissipation at the boundary leads to a sta- 
tionary state where the average weight saturates. In the 
language of interface dynamics this situation correspond 
with an elastic interface pinned at the boundary under 
an external force wq- In this case is better to look for a 
solution of the form 

w(x,t) = W(x) + y(x,t), (13) 

where W(x) is the solution of the stationary problem 

TV 2 W{x) + w = 0, (14) 

with homogeneous boundary conditions (W|boundary = 
0), and y(x,t) satisfies eq. ( pT| ) but with homogeneous 
boundary conditions. 

Again, y(x*,t) has zero average and, therefore, the 
mean weight in the stationary state is W(x). The so- 
lution of eq. ( |l4| ) can be easily obtained for certain ge- 
ometries. In d = 1 we can look for the solution in the 
interval < x < L with 1^(0) = W(L) = 0, obtaining 



W(x) 



w L 2 
2T 



x x 2 \ 
L^L 2 



(15) 



For d = 2, d = 3 and d > 3 (only the case d = 2 has 
a physical realization) we can look for the solution in 
a circle, sphere and hyper-sphere of radius L such that 
W = W(r), with < r < L and the boundary condition 
W(L) — 0, where r is the distance to the center. The 
solution is in this case given by 



W{r) 



w L 2 
2dT 



1 - 



L' 2 



(16) 



Thus, in any dimension, L is the characteristic length of 
the system and the saturation weight scales as L 2 . 

Now we are going to analyze the fluctuations around 
the average, described by eq. (11) with the correspond- 
ing boundary conditions and noise correlator. A formal 
solution of eq. ( p|) is given by 

y(x, t) = J df J d d x'G d (x, x>, t - t'; L)ri(x", t'), (17) 

where Gd(x, x', t; L) is the Green function of correspond- 
ing boundary problem, which depends on the dimension 
of the horizontal space d and on the linear size of the 
system L. The precise form of the Green function can be 
obtained only for some suitable geometries, however, in 
general it satisfies the scaling relation 



G d {x, x',l;L) 



1 



x x t 
L'T'L 2 



(18) 
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FIG. 1. Global average weight as a function of depth ob- 
tained from numerical simulations in one dimension. The 
continuous (dashed) line corresponds with absorbing (peri- 
odic) boundary conditions. In the case of absorbing boundary 
conditions the average weight saturates after a characteristic 
depth t s while it increases linearly with depth for periodic 
boundary conditions. 

where gd{x,x',t) depends on the spatial dimension and 
boundary conditions and L is a characteristic linear di- 
mension of the system, g^ can be obtained exactly, for 
instance, in one dimension and in a systems with radial 
symmetry. 

Using equations (jl7|), (|l8l), and (@) we compute the 
"surface roughness" Aw, obtaining the scaling relation 
in eq. ([)]) with z — 2 and a roughness exponent £ which 
depends on the choice of A(t). If the noise is uncorrelated 
in the vertical direction t (t c — > 0) then 



c 



2-d 



(19) 



which corresponds with the EW universality class. On 
the contrary, if it is strongly correlated in the vertical 
direction (t c — > oo) then 



c 



4-d 



(20) 



In both cases the exponent C, is independent of the choice 
of boundary conditions, and they are identical to the ex- 
ponents obtained from the Fourier analysis for a L — > oo 
system. Thus, while the average profile is strongly depen- 
dent on the boundary conditions, the "surface roughness" 
only depends on the noise correlator. 



III. NUMERICAL SIMULATIONS AND 
DISCUSSION 

To test results obtained in the previous section we have 
performed numerical simulations of the discrete model in 
one dimension (N = 2) with a uniform q- distribution. 
We have used w = 1 and lattice sizes L — 50, 100, 200 
and 400. Average were taken over 1000 realizations for 
the three smallest lattice sizes and over 100 realizations 
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FIG. 2. Weight profile in the stationary state for absorb- 
ing boundary conditions. The choice of scaled variables is 
explained in the text. The heavy line is a fit to the quadratic 
dependency P(x) — a + bx — cx 2 , with a — 0.02 ± 0.01, 
b = 3.91 ± 0.1 and c = 3.92 ± 0.1. 

for the largest one. Different magnitudes associated with 
force fluctuations were computed. 

Our analytical approach reveals that the existence of 
absorbing boundaries is necessary to obtain a satura- 
tion weight, otherwise the weight increases linearly with 
depth. In fig. [l] we have plotted the global average weight 



W(t) 



(21) 



at depth t, for periodic and absorbing boundary condi- 
tions. The agreement with our prediction becomes evi- 
dent. 

Moreover, in the case of absorbing boundary condi- 
tions, we have computed the average weight profile after 
saturation which in one dimension is given by eq. (|lE|). 
The average weight after saturation as a function of lat- 
tice position x is shown in fig. ^, for different lattice sizes. 
The scaled variables w/L 2 vs. x/L have been used as 
suggested by eq. (|l^) . All the curves collapse in a single 
plot showing that our choice of scaled variables is correct. 
Moreover, the scaled plot was fitted to the quadratic de- 
pendency P{x) = a + bx — cx 2 , with a = 0.02 ± 0.01, 
b = 3.91 ± 0.1 and c = 3.92 ± 0.1. These fitting parame- 
ters are in very good agreement with eq. (|l5|) since they 
satisfy a <C b and b ~ c. We have thus obtained the 
correct profile for the saturation weight. 

To investigate the transient region before saturation we 
have computed the global average weight as a function of 
depth. The result for different lattice sizes is shown in fig. 
||. The scaled variables W{t)/L 2 vs.t/L 2 have been used. 
For small depths W(t) scales linearly with t, as for peri- 
odic boundary conditions. In this region the force chains 
starting from bulk sites have not reached the boundary 
and, therefore, the weight dissipation at the boundary 
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FIG. 3. Average weight as a function of depth for different 
lattice sizes. The choice of scaled variables is explained in the 
text. 

is negligible. This transient behavior was previously ob- 
served in numerical simulations by Peralta-Fabi et al in 
one dimension Jt|. 

The average squared fluctuations (interface roughness) 
AW 2 = {(W — (W)) 2 ) as a function of t are shown in fig. 
g. Using the scaled variables AW 2 /L 2C , with C = 3/2, 
vs. t/L 2 we obtain a good data collapse. This scaled plot 
is very sensitive to the choice of £, a good data collapse is 
only obtained for £ = 1.50 ± 0.05. If the noise is uncorre- 
cted in the vertical direction then £ = 1/2, which is much 
smaller than our numerical estimate. On the contrary, if 
the noise is strongly correlated in the vertical direction 
then £ = 3/2, in very good agreement with our numerical 
estimate. When we avoid the multiplicative character of 
the noise we guess that, in compensation, strong correla- 
tions in the vertical direction should be considered. This 
supposition is now supported by numerical simulations. 

The simplicity of the continuum model we have pro- 
posed in the previous section allowed us to obtain some 
analytical results. We have avoided the complicated form 
of the noise, associated with the fluctuations in the frac- 
tions Qij, introducing an "annealed" noise with correla- 
tions in the vertical direction. Within this approxima- 
tion we have obtained the average weight in the station- 
ary state and characterized the fluctuations around the 
average, for both periodic and absorbing boundary con- 
ditions. 

In earlier numerical simulations by Liu et al ||, they 
considered periodic boundary conditions and compute 
the distribution of normalized weight v = w/t. They ob- 
serve that the distribution of normalized weight P(v) was 
independent of t for large t. This numerical result was 
later corroborated by Coppersmith et al using a MF the- 
ory for the discrete model |6[ . According to our approach, 
for periodic boundary conditions the average weight is 
given by wot/X and, therefore, (v) = wq/\. Since the av- 
erage normalized weight is independent of t so should be 
its distribution, in agreement with the numerical sim— 
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FIG. 4. Average squared fluctuations as a function of depth 
for different lattice sizes. The choice of scaled variables is ex- 
plained in the text. 

ulations and MF theory of the discrete model. 

The robustness of the asymptotic exponential decay 
of the normalized weight distribution P(v) is one of the 
main features of the g-model. We thus check if this be- 
havior is still observed when absorbing boundary con- 
ditions are considered. Before going to the numerical 
results we should remember that in this case for t ^> t s 
the average weight does not increases linearly with t but 
scales as L 2 . Hence, the appropriate choice of normalized 
weight is in this case v = w/L 2 . 

With this remark, in fig. (|^) we show the distribu- 
tion of normalized weight for different lattice sizes. First 
we notice that P(v) is independent of lattice size, sup- 
porting our choice of normalized weight. Moreover, the 
resulting universal curve clearly show an exponential de- 
cay for large v. However, for small v the data does not 
display the linear behavior predicted by the g-model for a 
uniform distribution . This discrepancy for small v is 
just a consequence of the different boundary conditions. 

On the other hand, in the case of absorbing boundary 
conditions our analytical approach is in agreement with 
previous MF theory by Peralta-Fabi et al |7) in one di- 
mension. Peralta-Fabi et al analyzed the case Qij = 0, 
obtaining the scaling dependencies W s ~ L 2 and t s ~ L 2 . 
We have shown that these scaling relations are valid in 
larger dimensions and are independent of the noise. How- 
ever, within the MF theory by Peralta-Fabi et al one 
cannot determine the scale dependency of the average 
squared fluctuations Aw 2 , which constitutes one of our 
main results. 

The scaling dependencies W s ~ L 2 and t s ~ L 2 are, 
nevertheless, in contradiction with the linear scaling ob- 
tained in classical theories § , some generalizations of the 
g-model 10 and even in experimental observations [ p^[ . 
However, this discrepancy cannot be attributed to our 
continuum analysis, which shows a very good agreement 
with the numerical simulations of the q-model. 
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IV. SUMMARY AND CONCLUSIONS 

We have obtained a coarse-grained equation of the dis- 
crete model introduced by Liu et al to describe force fluc- 
tuations in granular media. The multiplicative nature of 
the noise has been considered assuming strong correla- 
tions in the vertical direction. The stationary and tran- 
sient behavior were obtained analytically for periodic and 
absorbing boundary conditions. 

In this way we have demonstrated that the existence 
of walls, modeled by the absorbing boundary conditions, 
are a necessary condition to obtain an stationary state 
were the average weight is independent of depth. We 
have also shown that the the scaling of the saturation 
weight W s ~ L 2 and depth t s ~ L 2 with lattice size are 
valid in any dimension and arc independent of the noise, 
generalizing previous MF calculations in one dimension. 

For the first time we have obtained the scale depen- 
dency of the average squared fluctuations. The compari- 
son of this scale dependence with the one obtained in the 
numerical simulations support our guess about the exis- 
tence of strong correlations of the noise in the vertic 
direction. 
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FIG. 5. Distribution of normalized weight v = w/L 2 , after 
saturation <>t s , for the q- model with absorbing boundary 
conditions. The straight line is a fit to an exponential decay. 
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